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Abstract 

The Fortran package QCD-Pegasus is presented. This program provides fast, flexible 
and accurate solutions of the evolution equations for unpolarized and polarized parton dis- 
tributions of hadrons in perturbative QCD. The evolution is performed using the symbolic 
moment-space solutions on a one-fits-all Mellin inversion contour. User options include 
the order of the evolution including the next-to-next-to-leading order in the unpolarized 
case, the type of the evolution including an emulation of brute-force solutions, the evo- 
lution with a fixed number of flavours or in the variable-?7,j scheme, and the evolution 
with a renormalization scale unequal to the factorization scale. The initial distributions 
are needed in a form facilitating the computation of the complex Mellin moments. 



Program Summary 



Title of program: QCD-Pegasus 
Version: 1.0 
Catalogue identifier: 

Program obtainable from: http://arxiv.org/archive/hep-ph and its mirror sites by 
downloading the source of hep-ph/ 0408244 

Distribution format: uuencoded compressed tar file 

E-mail: avogt@nikhef.nl 

License: GNU Public License 

Computers: all 

Operating systems: all 

Program language: Fortran 77 

Memory required to execute: negligible (< 1 MB) 

Other programs called: none 

External files needed: none 

Number of bytes in distributed program, including test data etc.: 240 578 

Keywords: unpolarized and longitudinally polarized parton distributions, Altarelli-Parisi 
evolution equations, Mellin-space solutions 

Nature of the physical problem: solution of the evolution equations for the unpolarized 
and polarized parton distributions of hadrons at leading order (LO), next-to-leading order 
and ncxt-to-ncxt-to-leading order of perturbativc QCD. Evolution performed either with 
a fixed number of effectively massless quark flavours or in the variable-n^ scheme. The 
calculation of observables from the parton distributions is not part of the present package. 

Method of solution: analytic solution in Mellin space (beyond LO in general by power- 
expansion around the lowest-order expansion) followed by a fast Mellin inversion to 
x-space using a fixed one-fits-all contour. 

Restrictions on complexity of the problem: The initial distributions for the evolution are 
required in a form facilitating an efficient calculation of their complex Mellin moments. 
The ratio of the renormalization and factorization scales //f/A* has to be a fixed number. 

Typical running time: one to ten seconds, on a PC with a 2.0 GHz Pentium-lV processor, 
for performing the evolution of 200 initial distributions to 500 {x, ji) points each. For 
more details see section 6. 
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1 Introduction 



Parton distributions form indispensable ingredients for analyses of hard processes with 
initial-state hadrons, investigated in fixed-target experiments and at colliders like HERA, 
RHIC, TEVATRON and the forthcoming LHC. The task of determining these distribu- 
tions can, in principle, be divided in two steps. The first is the determination of the 
non-perturbativc initial distributions at some (usually rather low) scale (Iq. The second 
is the perturbative calculation of their scale dependence (evolution) to obtain the results 
at the hard scales fi. For the foreseeable future, the non-perturbative input cannot be 
calculated from first principles with a sufficient accuracy. Instead the initial distributions 
have to be fitted - and re-fitted once new data become available - using a suitable set of 
hard-scattering observables. In practice the evolution thus enters in both steps. 

Two methods for solving the evolution equations have been most widely applied in 
parton-distribution analyses. The first is the direct numerical integration of these integro- 
differential equations in {x, //.)-space, where x stands for the momentum fraction carried 
by the partons. Publicly available programs of this type can, for example, be found in 
refs. [1, 2] and [3]. In the second approach a MeUin transformation is apphed to turn 
the evolution equations into systems of ordinary differential equations (depending on the 
Mellin variable N) which are more easily accessible to a further analytic treatment. A 
C++ code of this type has been published in rcf. [4]. More programs of both types have 
been used and/or informally circulated in the perturbative-QCD community. 

In this article we present another A^-space evolution package. For easier reference 
the program has been given a name, QCD-PegasuS or PEGASUS in short, standing for 
'Parton Evolution Generated Applying Symbolic [/-matrix Solutions'. Here '[/-matrix' 
is a usual name for the key technical ingredient in the solutions which dates back, at 
least, to ref. [5]. The present program is a descendant of the code written by the present 
author fifteen years ago for the GRV analyses started in ref. [6]. Since then it has been 
almost completely rewritten more than once. Various intermediate versions have been 
regularly used in QCD studies by the author and by others. The present last incarnation 
of the program includes, in a now hopefully suflicicntly well documented and user-friendly 
manner, the evolution of unpolarized and helicity-dependcnt parton densities up to (in the 
former case) the next-to-next-to-leading order (N^LO) of perturbative QCD for any fixed 
ratio of the factorization and renormalization scales. The user can choose between various 
ways to truncate contributions of higher order (including emulations of the brute-force 
solutions) in both the fixed flavour-number and the variable-n^ evolution schemes. 

This manual is organized as follows. In section 2 we recall the formalism used for the 
A"-space evolution of the parton distributions. Topics related to the inverse MeUin trans- 
formations of the solutions and the initial distributions are then discussed in section 3. 
A compact user guide for the program is provided in section 4, followed in section 5 by 
a short reference manual of all routines. In section 6 we briefly address the accuracy and 
speed of the evolution by QCD-Pegasus, before we conclude in section 7. 
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2 The evolution equations and their solution 



In this section we discuss, in some detail, the formahsm employed in the program for the 
Mellin-space solution of the evolution equations. In the course of the discussion we point 
out some default choices made in the present version of the program, explain the main 
options available to the user, and indicate some structural restrictions. 



2.1 The running coupling constant 

As the reader will see below, the strong coupling constant plays a more central role in 
the present approach to the evolution of parton densities than usually in x-space programs. 
We therefore start the discussion with for which we employ the normalization 

At N™LO the scale dependence of ag is given by 

da 



dln//^ 



= /3N-Lo(as) = -5^ ag^+'/3fe , (2.2) 

k=0 



where Hr denotes the renormalization scale and stands for the number of effectively 
massless quark flavours, rij- is considered a fixed number until section 2.7. The expansion 
coefficients (3k of the /3-function of QCD are known up to /c = 3, i.e., N^LO 

/3o = 11 - 2/3 

(3i = 102 -38/371^ 

P2 = 2857/2 - 5033/18 ny+ 325/54 n/ 

Ps = 29243.0 - 6946.30 + 405.089 n/ + 1093/729 n/ . (2.3) 

Here the scheme-dependent quantities P2 [7, 8] and P3 [9] refer to the usual MS scheme. 
Only this scheme has been implemented in the program so far. For brevity the irrational 
coefficients of (3^ in Eq. (2.3) have been truncated to a sufficient accuracy of six digits. 

Eq. (2.2) can be integrated in a closed form only at low orders, and even then one only 
arrives at an implicit equation for as(/i^) beyond LO. The exact solution at NLO which 
expresses ag(/x^) in terms of its value ag(//o) at a reference scale /Iq, for example, reads 

^ +^.nf4VM.J"-!i;'!+!-°-S1!!l (2.4) 



Qsiiij) as{i4) V/^o/ las(/^o) [l + ^'ias(/^r)] 

with bk = Pk/Po- The program uses Eq. (2.4) only at LO (61 = 0), as in general a 
numerical iteration is required otherwise anyway. Beyond LO the value of as(A*r) by 
default determined directly from Eq. (2.2) by a fourth order Runge-Kutta integration [10]. 

All known orders m < 3 in Eq. (2.2) are available in the routine for Og to which m is 
transferred as NAORD. In the context of the evolution program, NAORD is not set externally 
but specified via the order chosen for the splitting functions, see section 2.2. 
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Another very common approach, used for instance by the Particle Data Group [11], 
is to expand the solution in inverse powers of La = ln(/i^/A^) where A is the QCD scale 
parameter. Up to N^LO this expansion yields [12] 



1 



1 



+ 



1 



bi In La 



1 



(PoLaY 
hi (-ln3LA + 



hi (in^ LA-lnLA-l)+62 



In^ La + 2 IuLa - - - 86162 IhLa + 



(2.5) 



Eq. (2.5) solves the evolution equation (2.2) only up to higher orders in I/La- As explained 
in section 2.4, this is an unwanted feature for the A^-space evolution which especially 
bedevils direct comparisons to x-space evolution programs. Therefore the use of Eq. (2.5) 
is not a standard option in the present evolution package. 



2.2 The general evolution equations 

The scale dependence of the parton distributions is governed by the evolution equations 

h{x, Pi') = P,j{x, Pi') fj{x, Pi') . (2.6) 

Here pi represents the factorization scale, and for the moment we put Hr = pi>- fi{x,pL') 
stands for the number distributions of quarks, antiquarks and gluons in a hadron, where x 
represents the fraction of the hadron's momentum carried by the parton. Summation over 
the parton species j is understood, and ® stands for the MeUin convolution. Eq. (2.6) 
thus represents a system of 2ny + l coupled integro-differential equations. The N™LO 
approximation for the splitting functions Pij{x,pb') reads 

m 
fc=0 

The splitting functions for the spin-averaged (unpolarized) case are now known at N^LO 
(= NNLO) [13, 14]. Note that Pi.j{x,fi') depends on pL only via the coupling as(yU^), a 
feature which forms the basis for the A^-space solution of Eq. (2.6) discussed in section 2.4. 

The splitting functions for the general case pir ^ pi can be obtained from Eq. (2.7) by 
Taylor-expanding as{pi') in terms of as(/i^). Up to N^LO this leads to 

Pij{x,i^,l^r) = as{nl)Pij\x) 

+ a'M){P^fi^)-PoPt,'k^)L) (2.8) 
+ iPS\^) - 2/?oLP,«(x) - {AL - (3',L'} pS\x)) 

+ a'M) (HPi^) - 3/3oLP,f (x) - {2AL - 3f3',L'} P^\x) 
- {/32L - 5/2 A/3oL2 + ^IL^] {x)) 
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with L = ln(/x^/ /x^). If L is a fixed number, then also in Eq. (2.8) the coefficients of Og (/x^) 
depend only on and the algorithms described below are applicable. In other words, the 
program is not designed to deal with choices like ji^ — + //^ where M is some mass 
scale. Note, however, that no such restriction is in place between physical scales and 

The order m in Eq. (2.7) (denoted by NPORD) and the ratio (denoted by FR2) are 

initialization parameters of the evolution package. The values m = 0, 1 and 2 are available 
at present for the standard MS factorization scheme. An extension to m = 3 — based 
on future partial results or even Fade estimates for P^^ — may be useful for uncertainty 
estimates in special cases, e.g., in determinations of Og from structure functions [15]. 



2.3 The flavour decomposition 

It is convenient to decompose the system (2.6) as far as possible from charge conjuga- 
tion and flavour symmetry constraints alone. The gluon-quark and quark-gluon splitting 
functions are flavour independent 

-fgq = -fgqi — -fgqi , -fqg = ''^f Pfiig, — "f^f Pqig ■ (2-9) 

Any difference Qi—qj and qi—qj of quark and (anti-) quark distributions therefore decouples 
from the gluon density g. Hence the combination maximally coupling to g is the flavour- 
singlet quark distribution 

Qs^i^i^r + qr) (2.10) 

r=l 



evolving according to 

dlnf,^ [g ) [P,q P^J'^yg ) ■ 
The singlet quark-quark splitting function Pqq is specifled in Eq. (2.15) below. 



(2.11) 



In order decouple the non-singlet (difference) combinations, we make use of the general 
structure of the (anti-) quark (anti-) quark splitting functions, 

PciiCik = Pmk = ^ikPqq + Pqq ■ (2-12) 

In general (beyond NLO), Eq. (2.12) leads to three independently evolving types of non- 
singlet combinations. The flavour asymmetries q^ and the total valence distribution q^^, 

(i^ik = Qi^Qi- {Qk ± qk) , qIs = Yl ((ir - qr) , (2.13) 

r=l 



respectively evolve with 

ns 

PZ = P^^-P^q + nf{P^^-P^^) ^ P^+P^ . (2.14) 



ns qq qq 
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Finally the singlet splitting function (2.15) can be expressed as 

Pqq = P„|" + nj(Pqq + Pqq) = P^^ + P^^ . (2-15) 

In the expansion in powers of flg, the flavour- diagonal ('valence') quantity Pqq in Eq. (2.12) 
starts at first order. Pqg and the flavour-independent ('sea') contributions Pqq and Pq^ — 
and hence the 'pure-singlet' term Ppg in Eq. (2.15) — are of order a^. A non-vanishing 
Pjfs ~ Pqq — Pqq iu Eq. (2.14) occurs for the first time at the third order. 

For the evolution of the fiavour asymmetries in Eq. (2.13) we use the basis 

k 

vf = ^{qi ± qi) - k{qk ± qk) (2.16) 

with k — 1, . . . , ny and the usual group-theoretical notation I = k"^ — 1. After performing 
the evolution, the individual quark and antiquark distributions can be recovered using 

+ = ^ - 7 + W^) - ('-'^^ 

where Vq = 0, together with the corresponding equation for the differences qi — qi- 



2.4 The AT-space solutions 

In the next two sections we describe the algorithm employed for the solution of the evolu- 
tion equations in Mellin-A'" space. Thus we now switch to the moments of all x-dependent 
quantities, 

a(N) = / dxx^-^a(x) . (2.18) 

The advantage of this transformation is that is turns the Mellin convolutions into simple 
products, 

[a^b]{N) = a{N)b{N) , (2.19) 

which greatly simplifies all further manipulations. The disadvantage of working in A^-space 
is that all quantities have to be known for complex values of N for the final transformation 
back to x-space. The resulting limitations of the program are discussed in section 3. 

As discussed above, we restrict ourselves to situations where the scale fi enters the 
right-hand side of Eq. (2.8) only through the (monotonous) coupling = as{fi^ = Kfi^) . 
Hence we can switch to ag as the independent variable. Using a matrix notation for the 
singlet system (2.11), the combination of Eqs. (2.2) and (2.7) yields 

— = {/3N-Lo(as)} PN--Lo{N,as) q{N,as) 



Poa. 
1 



+ al {P('\N) - hP^'\N) + {hi - h2)P^'\N)) + . . .] q{N, a^) 

oo 



k=\ 



q{N,a,) . (2.20) 



6 



In the last line we have introduced the recursive abbreviations 



R, = Lpi^), = ^P^^) -j^URu-i (2.21) 

for the splitting function combinations entering the new expansion (2.20). hk has been 
defined below Eq. (2.4), and P^^'^ is the coefficient of 0^^+^ in Eq. (2.8). As in Eq. (2.21), 
we will often suppress the explicit reference to the Mellin variable N below. 

The singlet sphtting- function matrices Rk of different orders k do not commute. There- 
fore the solution of Eq. (2.20) cannot be written in a closed exponential form beyond LO. 
Instead we employ a series expansion around the lowest order solution. 



q^o{N, as, iV) = M q(iV, oq) = L{N, a,, oq) q{N, oq) (2.22) 

Van/ 



-Ro{N) 

with oo = as(/i^o = /t/io)- This expansion reads 

q{N,a,) = UiN,a,)L{N,a,,ao)U-\N,ao)q{N,ao) (2.23) 



k=l 



l + X]«st7fc(iV) L{a,,ao,N) l + J^^UkiN) 



k=l 



q{ao,N) . 



Here the third, a^-independent factor normalises the evolution operator to the unit matrix 
at /Xq, instead of to the LO result (2.22) at infinitely high scales. The evolution matrices 
Uk are constructed from the splitting function combinations Ri<k in the next section. 

Eqs. (2.20) and (2.23) offer various ways to define the N™LO solution which differ in 
terms of order n > m. In the program the choice between the resulting options is made by 
the initialization parameter IMODEV. One obvious choice is to keep the terms originating 
from /^N-^LO in Eq. (2.2) and Pn^lo in Eq. (2.7) to all orders (in practice: to a sufficiently 
high order) in both (2.20) and (2.23). This is equivalent to a direct iterative solution of 
Eq. (2.6) as performed by standard x-space evolution programs, to which the results can 
then be compared directly. This mode for the evolution is invoked for IMODEV = 1. 

Note that this equivalence only holds if as(/x^) solves Eq. (2.2) exactly. Otherwise, for 
example if Eq. (2.5) is used, the evolution equation (2.20) in Og is obtained from Eq. (2.6) 
by dividing the l.h.s. and the r.h.s. by quantities which differ somewhat, viz das/dln /j!^ 
and /SNmLO- The difference introduced by this mismatch is a higher-order effect, but its 
numerical impact is far from negligible, e.g., most of the difference between the two parts 
of table 1 in the 1996 NLO comparisons [16] arises from this source. 

One can equally well take the point of view that at N"'LO terms beyond a™ should be 
removed in the square brackets in Eq. (2.20), as these terms would receive contributions 
from p(">™) and Pn>m- At N^LO, for example, one then retains only the terms explicitly 
written down in the second and third fine of Eq. (2.20). If still 'all' orders are kept in 
the solution (2.23), one arrives at a second iterative option employed for IMODEV = 2. 
Finally one can instead apply the same reasoning to the matrices Un>m in Eq. (2.23) 
which would also receive contributions from and Pn>m- Expanding also the 
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term in Eq. (2.23), necessary in the singlet case as explained below Eq. (2.31), one then 
arrives at the so-called truncated solution. Up to N^LO this solution is thus given by [17] 

QnsloM = L + asUiL - aoLUi 

+ a^^U2L-a,aoUiLUi + alL[u^ - U^) 

+ al UsL- alao U2LU1 + a,al Ui L [uj - U2) 

-alL(uf- U,U2- U1U2+ t/3)] q(ao) , (2.24) 

where we have suppressed all arguments of L{N, Cg, Oq) for brevity. The NLO (NNLO) 
approximations are obtained from Eq. (2.24) by respectively retaining only the first (first 
and second) line in the square bracket. These truncated N™LO solutions (at present 
implemented for m < 2) are employed by the program for any value of IMODEV other than 
1, 2 and, for the non-singlet cases, 3. The latter case is addressed in section 2.6. 

The approaches discussed above obviously differ only in terms beyond the order under 
consideration. The iterative procedures introduce more scheme-dependent higher-order 
terms into the evolution of observables in a general factorization scheme. On the other 
hand, the truncated solution does not satisfy the evolution equations (2.6) literally, but 
only in the sense of a power expansion, i.e., up to terms of order n > m like Eq. (2.5) for 
the coupling constant. The differences between the respective results can be regarded as 
a lower limit for the uncertainties due to higher-order terms. 



2.5 The evolution matrices Uk 



Inserting the ansatz (2.23) into the evolution equations (2.20) and sorting in powers of 
tts, one arrives at a chain of commutation relations for the expansion coefficients Uk{N) : 



[U^,Ro\ = Rl+U^ 

[U2,Ro] = R2 + R1U1 + 2U2 



(2.25) 



fc-i 



[Uk,Ro] = Rk + J2^k-iUi + kUk = Rk + kUk 



i=l 



For the fiavour-singlet system (2.11) these equations can be solved recursively by applying 
the eigenvalue decomposition of the LO splitting function matrix [17, 18], completely 
analogous to the classical truncated NLO solution with only Ui in ref . [5] . One writes 



Rq = r_e_ + r+e+ , 
where r_ (r+) stands for the smaller (larger) eigenvalue of Rq, 

1 



2/3o 



(0) 

gg 



(2.26) 



(2.27) 
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The matrices e± denote the corresponding projectors, 

1 

r± — Tip 

and / the 2x2 unit matrix. Thus the LO evolution operator (2.22) can be represented as 



(2.28) 



L{as,ao,N)^e4N)( 



.(JV) 



+ e. 



-iN)( 



-r+(iV) 



(2.29) 



Inserting the identity 

Uh = e_ C4e_ + e_ £7^6+ + e+ C4e_ + e+ C4e+ (2.30) 
into the commutation relations (2.25), one finally obtains the coefficients in Eq. (2.23), 



1 r ~ ^1 e.RkB 
Uk = -7 e_iifee_ + e+Rke+ H 



k r I 



k 



(2.31) 



Note that the poles in Uk{N) at iV- values where r_(A^) — r+(A^)±A; vanishes arc cancelled 
by the U^^ term in the solution (2.23). The expansion of U^^ mentioned above Eq. (2.24) 
achieves this cancellation also for the truncated solutions. This can be made directly 
visible by inserting the explicit forms (2.29) and (2.31) into the solution (2.24) and writing 
the contributions in an appropriate order. At NLO, for example, one arrives at 



e_ + (ao — Os) e^RiB. 



(2.32) 



ao — Qs 



aoJ 



-r+ 



e-Rke+ 
r+ — r_ — 1 



+ ( + ^ - ) ^ q(«o) • 



Not only the denominator in the second line vanishes for r+ — r_ = 1, but also its 
coupling-constant prefactor in the round brackets. 



2.6 Non-singlet cases and symmetry breaking 



Eq. (2.25) also holds for the scalar evolution of the non-singlet combinations (2.13) of the 

quarks distributions, but with the obvious simplification that the right-hand sides vanish. 
This facilitates a direct recursive solution for in which, unlike the singlet results (2.31), 
no spurious poles occur. Consequently the truncated solution can be written down also 
without the expansion of in this case, at N™LO yielding 



^ns 



(«s) 



k=l 



-±,v 



k=l 



^ns 



(ao) 



(2.33) 



This solution is accessed by IMODEV = 3 (together with Eq. (2.24) for the singlet case). 



Both iterated non-singlet solutions can be written down in a compact closed form at 
NLO. Hence instead of the expansion (2.23) we use for IMODEV = 1 

bi \1 + biaoj J Voo. 



?NLo(as) = exp 



do) , 



(2.34) 
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and for IMODEV = 2 

_ ^ ns 

?NLo(«s) = exp[{a,-ao)U^] I^^J ° g„=^(ao) . (2.35) 

At NNLO the non-singlet solutions have been programmed analogous to the singlet case. 

Due to the differences of P^t -fnl for m > 1 and of and for m > 2 in 
Eq. (2.14), qualitatively new effects arise at the respective orders m, namely a breaking 
of symmetries imposed on the initial distributions. 

Consider the NLO evolution of an input u = + u, d = dy + d with an SU(2)- 
symmetric sea, m(/Uq) = d{^'^). The initial distributions for the evolution of and ^3" in 
Eq. (2.16) are then identical, but due to ^ Ri not the results of the evolution, 

{uy - (iv)(as) 
v+(as) 

Subtracting these two (truncated) solutions yields 

2(u-d){a,) = (a,-ao) - (^j ° {u, - dy){ao) . (2.37) 

Hence a flavour symmetry of the input sea quark densities is not preserved by the NLO 
evolution, if the valence distributions (as in the case of the proton) are not identical. The 
same procedure applied to and in Eq. (2.13) shows that s 7^ s at NNLO even for 
(s — s)(/io) = 0. Both effects are very small. The reader is referred to ref. [19] for a further 
discussion of especially the strange-sea asymmetry. 

2.7 The treatment of heavy flavours 

The evolution of the strong coupling and the parton distributions can be performed in 
both the fixed flavour-number scheme (FFNS) and the variable flavour-number scheme 
(VFNS). This choice is made via the initialization parameter IVFNS. The former scheme 
is used for IVFNS = 0, any other value leads to the latter. The number of flavours Uj- for 
the FFNS evolution is specified by the initialization parameter NFF. The values NFF = 
3 ... 6 can be used. For IVFNS 7^ the number assigned to NFF is irrelevant. 

In the VFNS case we need the matching conditions between the effective theories with 
rij and Uj + 1 light flavours for both the strong coupling and the parton distributions. 
For as these conditions have been determined at N^LO and N^LO in refs. [20] and [12], 
for the unpolarized parton densities they are known to N^LO from ref. [21]. 

In the present program the transitions + 1 are made, for both ag and the 

parton densities, when the factorization scale equals the pole masses of the heavy quarks, 

^ ml , h ^ c, b,t . (2.38) 



{l + (a,-ao)i?r}f- 



= + (ds - ao) Rf] ( — ] {uy-dy){ao) . 



dy){ao) 



\aoJ 



(2.36) 
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For this choice the matching conditions for the parton distributions read up to N™ ^LO 
lt^^'\x, ml) = lt^\x, ml) + 5m2 A;if\x) ® li''^\x, ml) (2.39) 
where I — q, q and i — 1, . . .rif, and 

^("/+^)(x,m2) = g(-^f\x,ml)+ (2.40) 
Sm2 al [Af^;? {x) ® (x, ml) + A'^f (x) ^("/^ (x, ml) 

{h + h)^^f+'\x,ml) = [i^f (a;)®gf^\a;,m2) + i^f (x)®^7("/)(a;,mD 

with h = h. The coefficients A*^^^ for the spin-averaged case can be found in Appendix B 
of ref. [21] from where also their notation has been taken over. Due to our choice (2.38) 
for the thresholds only the scale-independent parts of the expressions for A^'^^ are needed. 



[KmD = as"^\Kml) + J2 ■ (2.41) 



The corresponding N™LO relation for the coupling constant at /x^ = k/j,^ is given by 

n=l 1=0 

The pole-mass coefficients c„<3^; in Eq. (2.41) can be inferred from Eq. (9) of Ref. [12], 

where 4ai ■'^ ^ is expressed in terms of 4ai ^\ i.e., the matching is written backward in 
rif for a different normalization of the coupling. In our notation these coefficients read 

2 

^1,0 ~ 5 '^1,1 ^ 



3 

_ 14 _ 38 _ 4 

"2,0 ~ ~^ 5 '^2,1 ~ ~^ ) — - 



(2.42) 



and 



C3 = 340.729 - 16.7981 

8941 _ 409 _ 511 _ 8 

^3,1 - 27 ' '^^'^ ~ ~9~ ' '^^'^ ~ 27 ■ ^ ^ ^ 

As in Eq. (2.3), we have truncated the irrational N^LO coefficients to six digits here. 
Note that, while the parton distributions are continuous at the flavour thresholds for our 
choice (2.38) up to NLO, the same holds for the NLO coupling constant only under the 
additional condition yU,. = /i leaving only the vanishing coefficient c^q in Eq. (2.41). At 

NNLO we use ai ■'^ \Kml) on the right-hand sides of Eqs. (2.39) and (2.40). 

If the program is run in the variable flavour-number mode, the initial distributions 
are specifled for = 3, and Eq. (2.40) is employed at least for the charm distributions. 
Therefore the initial factorization scale /iq has to satisfy fiQ ^ n^-c and, consequently, not 
too small a value should be chosen for k. The inclusion of top (and bottom) among the 
partons can be switched off by simply assigning a sufficiently large value to the respective 
mass. The masses m^, like /iq, the initial coupling as(//o) the initial light-parton 
distributions, are not fixed at the initialization of the program, but are defined (and can, 
of course, be re-defined) at a later stage as explained below. 
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3 Complex moments and the Mellin inversion 



In this section we discuss the issues related to the inverse Melhn transformation required 
for recovering the x-spacc parton distributions (and, in general, related observables) from 
the A'^-space expressions used for the intermediate calculations. The discussion includes 
our choice of the Mellin-inversion contour, the required analytic continuations and the 
restrictions on the form of the initial distributions resulting from our approach. 

3.1 From moments to ic-space 

We first consider the inverse transformation of the Mellin moments (2.18). If, as in our 
cases, a{x) is piecewise smooth for x > 0, the corresponding Mellin inversion reads 

1 rc+ioo 

a{x) = — / dNx-^a{N) , (3.1) 

where the real number c has to be such that dx x'^~^a{x) is absolutely convergent [22]. 
Hence c has to lie to the right of the rightmost singularity iVmax of a{N). The contour 
of the integration in Eq. (3.1) is displayed in Fig. 1 and denoted by Cq. Also shown is 
a deformed route Ci, yielding the same result as long as no singularities Nj, of a{N) are 
enclosed by Cq — Ci. The Ni are real with Ni < A^max < c for the N™LO evolution of 
parton distributions, thus this requirement is fulfilled automatically in our case. 
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Figure 1: Two integration contours for the inverse Mellin transformation in Eq. (3.1). 
The crosses schematically indicate the singularities of the moment-space function a{N). 

It is useful to rewrite Eq. (3.1) as an integration over a real variable. We are deahng 
with functions which satisfy a*{N) — a{N*), where * denotes the complex conjugation. 
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In this case contours characterized by the abscissa c and the angle as in Fig. 1 yield 

xaix) ^ - d^Im e^'^a;'-'=-"""P(^'^)a(Ar = c + ^e^'^) . (3.2) 

While formally the integral does not depend on c > A^max and 7r>0>7r/2,a suitable 
choice of these parameters is obviously useful for an efficient numerical evaluation. In 
particular, it is advantageous to choose > 7r/2 since then the factor exp ^zlog^cos^^ 
introduces an exponential dampening of the integrand with increasing in addition to 
the oscillations already present for = 7r/2, i.e., the textbook contour Cq. Consequently 
a smaller upper limit Zraax can be employed in the numerical implementation of Eq. (3.2). 

The contour implemented in the present program forms a practitioner's compromise 
between speed and accuracy for the cases encountered in QCD evolutions. The same 
contour is employed for all inversions, and on that contour a fixed set of intervals is 
defined. For each of these intervals we perform a standard eight-point Gauss-Legendre 
integration [10]. As said, all this is done regardless of x, the initial distributions and all 
evolution parameters discussed in the previous section. Hence all moments of the splitting 
functions, evolution-operator coefficients and matching conditions need to be calculated 
only once on the resulting fixed grid of complex A^-values at the initialization of the 
program. Note that this approach is rather different from that of refs. [23, 4], where a 
parabolic contour is optimized with respect to the shape of the initial distributions. 

Specifically, we choose (j) = 3/4 tt (too close to neither Cq nor the singularities of the 
integrands) and c = 1.9 (see below). The upper limit ^max depends on x in accordance 
with the ln(l/a;) damping mentioned below Eq. (3.2) above: Eight intervals covering the 
region < z < 5 are used for x < 0.01 (of which four have < z < 1). Three more 
intervals with 5 < z < 14 are added for 0.01 < x < 0.3, and another three covering 
14 < z < 32 are included for 0.3 < x < 0.7. Above the latter value of x the final four 
intervals with 32 < 2; < 80 are included as well. Under these conditions the chosen value 
c = 1.9 for the abscissa is a compromise between high accuracy for very small quantities at 
very large x (improved for larger values) and very small x (improved for smaller values). 

This standard setup, used for the default initialization parameter IFAST = is suffi- 
cient for a five-digit accuracy of the evolution of the proton structure at 10^^ ^ x < 0.9, 
with the exception of the tiny antiquark distributions at a; ~ 0.9. A yet faster, but at very 
large and very small a;- values less reliable inversion can be invoked by IFAST 7^ 0. The 
program then runs with minimally 4 and maximally 10 (instead of 8 and 18) intervals. 



3.2 Splitting functions and matching coefficients 

The complex- iV moments required for the evolution cannot be computed by direct numer- 
ical integrations of Eq. (2.18) on our tilted contour Ci shown in Fig. 1. Therefore we have 
to work with the proper analytic continuations or, where these are not known, with suffi- 
ciently accurate x-space parametrizations of which the moments are known analytically. 
That situation already occurs at NLO, but becomes more severe at higher orders. 
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The unpolarized NLO splitting functions have been programmed, long ago, using the 
integer- A?" results in Eq. (2.5) of ref. [24] (see also section 5 of ref. [25]) and in Appendix 
B of ref. [26], together with the analytic continuations provided in Appendix A.l of 
ref. [6]. With one exception, these analytic continuations can be expressed in terms of 
the complex polygamma functions ^^"-\z). These functions are calculated from their 
asymptotic expansions [27], after for |Im(z)| < 10 using the functional equations to shift 
the argument to Re(^) > 10. The exception are the moments 

f'dxx--'^ (3.3) 
Jo 1 + x ^ ^ 

where Li2 is the standard dilogarithm. Improving on Eq. (A. 6) of ref. [6], these moments 
can be approximated to a sufficient accuracy by inserting the parametrization 

^ 1 - 0.9992 x + 0.9851x2 -0.9005^3 

1 + x 

+ 0.6621 - 0.3174 x^ + 0.0699 (3.4) 

into Eq. (3.3) and using the standard expression [28] for the moments of Li2{x). A useful 
hst of Mellin transforms can also be found in the appendix of ref. [29] , see also ref. [30] . 
The polarized NLO splitting functions [31, 32, 33] involve the same set of functions; they 
have been implemented using the appendix of ref. [34] . 

The recent complete integer- expressions for the unpolarized NNLO sphtting func- 
tions [13, 14] in moment space include harmonic sums [35] up to weight five, of which 
the analytic continuations are not known at present. We therefore resort to the accurate 
x-space parametrizations of these functions provided in Eqs. (4.22) - (4.24) of ref. [13] 
and Eqs. (4.32) - (4.35) of ref. [14]. The complex moments of these approximations 
can be expressed in terms of the polygamma functions. Based on the accuracy of these 
parametrizations and the size of the NNLO effects in the evolution [13, 14], we expect 
that the relative errors introduced by this procedure amount to about 10~^ or less. 

Finally we need, for the NNLO variable flavour- number evolution, the complex mo- 
ments of the flavour-matching coefficients A^"^^ of ref. [21]. These can be expressed in 
terms of the psi-functions (-2) as well, with the exception of A^^'^\x) in Eq. (2.40). 
Also for this function the program presently uses the moments of a parametrization, viz 

A^^^\x) ^ - 1.111 ln-'^(l - x) - 0.400 ln2(l - x) - 2.770 ln(l - x) 

- 187.8 + 249.6 x - 146.8 In^ x ln(l - x) - 93.68 \nx 

- 3.292 W X - 1.556 In^ x - 24.89 x'^ - 0.006 S{1 - x) . (3.5) 

The convolution of the approximation (3.5) with typical gluon distributions differs from 
the exact results by less than one part in a thousand except close to zeros of these results. 
Note that the parametrization (3.5) has been employed for the approximate [36] VFNS 
NNLO part of the 2001/2 evolution benchmarks presented in table 6 of ref. [37]. 



14 



3.3 The initial distributions 



Also the initial conditions fi{x, /Xq) for the evolution are, of course, needed by the program 
in a form which facilitates a computation of the moments on the Mellin inversion contour. 
This is, presumably, the most severe restriction of the flexibility of the A^^-space approach 
as compared to direct numerical x-space solutions. For the latter, all quantities are usually 
defined on discrete grids of x-values, hence no functional form whatsoever is required for 
fi{x, /^o)- Foi' the present program a functional form is definitely needed, and furthermore 
that form should preferably be such that its complex moments can be readily computed. 
This is, for example, not the case for the form employed in the CTEQ 6 analysis [38] , 



'1 + AixY' e 



dtx 



(3.6) 



An A^-space evolution of Eq. (3.6) would require an accurate internal re-par ametrization. 

What can be readily handled by an A^-space program, on the other hand, should be 
perfectly adequate for a sufficiently general ansatz for the initial distributions. In the 
present version of the program, we have included the two six-parameter standard forms 



i,5 



xPiA Jf- p.^^ X 



and 



xfiix,IJ^t) = Ar,p.,x^i,2(i_a;)Pi,3 I i+p.^a;0.5^^.^^^^.^ 



1.5 



(3.7) 



(3.8) 



The moments of these functions are given in terms of Euler's Beta function B{zi,Z2). 
This functions is implemented using the asymptotic expansion of the logarithm of the 
Gamma function (the Stirhng formula) [27] at Re(2;i) > 5 and Re(2;2) > 5 together with 
the functional equation. The A^-dependent argument zi is first inverted for Re(2;i) < —10, 
i.e., the asymptotic expansion is invoked for zl = 1 — zi — Z2- 

Speed is a much more important issue here than for the initialization of the splitting 
functions, [/-matrix elements and matching conditions: In fits of the parton densities, we 
need to deal efficiently with a large number of calls of the initial distributions, each of 
which in turn requires on the order of 10^ evaluations of B{zi,Z2). Note that the last 
term in the square bracket in Eq. (3.7), like the corresponding last two terms in Eq. (3.8), 
does not require new calls of B{zi, Z2) since the functional equation in the first argument 
can be used instead. Therefore the extension of Eqs. (3.7) and (3.8) to (many) suitably 
chosen higher powers of x does not pose any efficiency problems. 

The ansatz (3.7) or (3.8) is used for the initial u and d valence-quark distributions 
— u — — d — the corresponding antiquark (sea) densities L+ = 2{u -\- d) 

and L_ — d — u, the gluon distribution g, and for the strange-fiavour combinations 

s± = s ± s. Heavy-flavour initial distributions h{x,fiQ) are obviously not required for the 
VFNS evolution starting with = 3, see section 2.7. At present, independent shapes 
for h{x, yUp) are not included for the FFNS evolution either. The definitions and available 
user-options for the extra coefficients Aj in Eqs. (3.7) and (3.8) will be explained below. 
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4 A brief Qcd-Pegasus user guide 



Most of the routines building up the evolution package are not directly called by the 
user for standard applications. All he/she needs to interact with, are the routines for the 
initialization of the program, the specification of the initial distributions and parameters 
like as(A*o)' reading out the results of the evolution. In this section we describe the 

available input and output variables of these routines and show a small example program. 

4.1 General initialization 

All input-independent quantities required for the unpolarized evolution (the polarized 
case is deferred to section 4.5) and n^-matching described the section 2 are initialized by 

CALL INITEVOL('EVOLPAR') . 

The integer parameter EVOLPAR defines how values are assigned to the six initialization 
parameters NPORD, FR2 introduced in section 2.2, IMODEV discussed in sections 2.4 and 
2.6, IVFNS and NFF of section 2.7 and IFAST defined in section 3.1. For EVOLPAR = 1, 
INITEVOL reads these parameters from a six-line file usrinit . dat looking like 

IFAST 

1 IVFNS 
4 NFF 

1 IMODEV 
1 NPORD 
1 . ODO FR2 . 

For EVOLPAR = 2, INITEVOL obtains the corresponding values by calling the subroutine 

USRINIT (IFAST, IVFNS, NFF, IMODEV, NPORD, FR2) 

provided by the file usrinit . f . For any other other value of EVOLPAR, the program 
uses a default set of initialization parameters, actually the values displayed above. 

The available initialization options can be briefly summarized as follows : 

IFAST 

Speed/accuracy flag for the Mellin inversion. IFAST = is the standard. Any other value 
leads to a faster, but especially for very large and very small x less reliable version. 

IVFNS 

Switches between the evolution with a fixed number of fiavours (for IVFNS = 0) and that 
in the variable flavour- number scheme, starting with = 3 (for any other value). 

NFF 

The fixed number rif of fiavours (between 3 and 5) for IVFNS = 0. Unused for IVFNS ^ 0. 
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IMODEV 

Switches between the various evolution modes beyond LO. IMODEV = 1 and 2 invoke the 
respective iterated solutions of the evolution equations truncated in d/dln/j,'^ and d/dttg. 
Other values (IMODEV = 3 is special) lead to the faster truncated solutions. IMODEV = 1 
emulates the standard x-space treatment if as(A*r) adequate, i.e., not truncated itself. 

NPORD 

The perturbative order m = 0, 1 or 2 of the evolution, defined as the 'm' in N°^LO. 
FR2 

The constant ratio /^^//i^ of the factorization scale /i and the renormalization scale /ir- 



4.2 Input parameters and initial distributions 

The input parameters and initial light-parton distributions for the evolution are set by 

CALL INITINP('INPPAR') . 

Analogous to EVOLPAR in the previous section, the integer parameter INPPAR specifies how 
the input parameters are read in. If INPPAR is neither 1 nor 2, the program will evolve a 
default input, viz the one used for the 2001/2 benchmark tables [37], 

xu^{x,iJ,l) = 5.107200a;°-^(l-x)^ 
xd^{x,f4) = 3.064320 (1 -x)^ 

xg{x,nl) = 1.700000 - (4.1) 

xd{x,f4) = . 1939875 ,T"°-^(1 -a;)'^ 

xu {x, nl) — (1 — x) xd{x, iiq) 

xs{x,iil) — xs{x,iil) — 0.2 x{u + d){x, III) 

with 

a,{iJ,l^2 GeV^) = 0.35 (4.2) 
and, in the variable flavour-number case 

rric ^ Ijio, mb = 4.5 GeV , rrit = 175 GeV . (4.3) 

INITINP reads these parameters from a file usrinp.dat for INPPAR = 1. For the above 
input (using the ansatz (3.7) for definiteness) that file may look like 

2 . ODO M20 

. 35D0 ALPHSI 

2 . ODO MC2 

20 . 25D0 MB2 

3 . 0625D4 MT2 

1 NFORM 
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1 



IMOMIN 
ISSIMP 
0.8D0, 
0.8D0, 
. 9D0 , 
-O.IDO, 
0.9D0, 
-O.IDO, 
-O.IDO, 







2.0D0, 
l.ODO, 



3.0D0, 
4.0D0, 
6.0D0, 
6.0D0, 
6.0D0, 
6.0D0, 
5.0D0, 



0.5D0, 
0.5D0, 
0.5D0, 
0.5D0, 
0.5D0, 
. 5D0 , 
0.5D0, 



O.ODO, 
O.ODO, 

O.ODO, 
O.ODO, 
O.ODO, 
O.ODO, 
O.ODO, 



O.ODO 
O.ODO 
O.ODO 
0.5D0 



PUV 
PDV 



0.193987D0, 
0.136565D0, 
. ODO , 
0.027313D0, 
l.ODO, 



O.ODO 
0.5D0 



O.ODO 



PLM 



PLP 



PSM 
PSP 



PGL 



For INPPAR = 2, INITINP obtains the corresponding values by calling the subroutine 
USRINP (PUV, PDV, PDL, PDL, PLS, PSS, PGL, M20, ALPHSI, 



provided, for example, by the file usrinp . f . This subroutine can also be modified to 
form a fit-parameter interface to programs like MiNUlT. 

The input parameters and options for the initial distributions presently available are : 

M20 

The initial factorization scale /ig- Equal to or smaller than for the VFNS evolution. 
ALPHSI 

The strong coupling [not our internal Og = Q!s/(47r)] at M20 = /^q (even for /ij. ^ /J,). 
MC2 < MB2 < MT2 

The squared c, b and t masses for the VFNS case IVFNS 7^ 0. Irrelevant for IVFNS = 0. 
NFORM 

Switches between the functional forms (3.7) for NFORM = 1, and (3.8) for NFORM = 2. 
IMOMIN 

For IMOMIN = 0, Nl_^ = Ng^ = 1 is used in Eqs. (3.7) and (3.8). Otherwise these factors are 
such that p^^^i = PLP(l) andp^^^^ = PSP(l) are the corresponding momentum fractions. 

ISSIMP 

Switches between the full input ansatz (for ISSIMP = 0) and a simplified input (otherwise), 
s + s = PSP(1)-L_|_ with s — s, for the initial strange-fiavour distributions. 

PUV, PDV, PLM, PLP, PSM, PSP, PGL 

The dimension-6 arrays for the input parameters p^j in Eqs. (3.7) and (3.8) for Uy, dy, 

L_ = d — u, = 2(d + ?i), s± = s±s and g. PUV(l) and PDV(l) are the respective quark 
numbers (first moments), 2. DO and l.DO for the case of the proton. PGL(l) is the total 
fractional momentum carried by the partons (usually l.DO). Nl_ = Ns_ = 1. PSM(5) is 
not used, the corresponding input parameter is fixed by the vanishing of the first moment. 



MC2, MB2, MT2, NFORM, IMOMIN, ISSIMP) 
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4.3 The aj-space parton distributions 

The evolved parton distributions Mellin-inverted to x-space are accessed by 

CALL XPARTON (PDFX, AS, 'X', 'M2^ 'IFLOW, 'IFHIGH', 'IPSTDO. 

Here PDFX is the double-precision parton output array, user-declared before as PDFX (-6 : 6) . 
AS (also double-precision, as all real variables) returns the corresponding value as(/i^(/i^)). 
X and M2 are simply the x and fi"^ for which the evolved partons are requested from 
XPARTON. The other three (integer) parameters are flags helping to avoid wasting time. 
For IPSTD = 0, the notation for the index of PDFX is 

PDFX(O) = 5(, PDF(l)=Mv, PDFX(-l) = u + u, PDF(2) = ... 

and otherwise 

PDFX(O) = g, PDF(l) = u, PDFX(-l) = u, PDF (2) = d, .... 

Note that XPARTON return xg etc. IFLOW and IFHIGH are lower and upper limits for the 
values of this index for which the Mellin-inversion is actually performed. For example, if 
the program is run for = 3 parton flavours, there is no point (except once, for checking) 
in letting the program work to produce the inevitable zeros for c{x,ii^) etc, thus IFLOW 
= —3 and IFHIGH = 3 should be used. If, moreover, the user is interested (this time) 
only in the combinations qi + (ji and the gluon density, then IPSTD = 0, IFLOW = —3 and 
IFHIGH = would be most efficient. Also the calculation of the valence distributions is, 
in terms of speed, better performed using IPSTD = 0. 

4.4 An example program with output 

We now present a small main program (provided by the file IhOltab.f ) which can be 
used to check whether the evolution package is working properly at least up to NLO. The 
internal standard inputs are used in the calls of both INITEVOL and INITINP. 

PROGRAM LHOITAB 

* 

IMPLICIT DOUBLE PRECISION (A - Z) 
DIMENSION PDFX(-6:6), XB(ll) 
PARAMETER ( PI = 3.1415 92653 58979 DO ) 
INTEGER Kl 

* 

* . .Access the input parton momentum fractions and normalizations 

COMMON / PANORM / AUV, ADV, ALS, ASS, AGL, 

NUV, NDV, NLS, NSS, NGL 

* 

* . .The values for mu"2 and x 
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DATA M2 / 1.D4 / 

DATA XB / l.D-7, l.D-6, l.D-5, l.D-4, l.D-3, l.D-2, 
l.D-1, 3.D-1, 5.D-1, 7.D-1, 9.D-1 / 

. .General initialization (internal default) 
CALL INITEVOL (0) 

..Input initialization (internal default) 
CALL INITINP (0) 

. .Output of the momentum fractions and normalizations 

10 FORMAT (2X,'AUV ,F9.6,2X, 'ADV =' ,F9 . 6 , 2X, ' ALS =',F9.6, 

2X,'ASS =' ,F9.6,2X, 'AGL =\F9.6) 

WRITE (6, 11) NUV, NDV, NLS, NSS, NGL 

11 FORMAT (2X,'NUV =' ,F9.6,2X, 'NDV =' ,F9.6,2X, 'NLS =',F9.6, 

2X,'NSS =' ,F9.6,2X, 'NGL =',F9.6,/) 

. .Loop only over x (M2 is fixed in the Les-Houches tables) 
DO 1 Kl = 1, 11 
X = XB (Kl) 

. .Call of the Mellin inversion, reconstruction of L+ and L- 
CALL XPARTON (PDFX, AS, X, M2, -5, 2, 0) 

LMI = (PDFX(-2) - PDFX(-l) - PDFX(2) + PDFX(l)) * 0.5 
LPL = PDFX(-l) + PDFX(-2) - PDFX(l) - PDFX(2) 

. .Output to be compared to the upper part of Table 4 
IF (Kl .EQ. 1) WRITE (6,12) AS * 4*PI 

12 FORMAT (2X, 'ALPHA_S(MR2(M2)) = ',F9.6,/) 

IF (Kl .EQ. 1) WRITE (6,13) 

13 FORMAT (2X, 'x' ,8X, 'xu_v' ,7X, 'xd_v' ,7X, 'xL_-' ,7X, 'xL_+' , 

7X, 'xs_+' ,7X, 'xc_+' ,7X, 'xb_+' ,7X, 'xg' ,/) 

WRITE (6,14) X, PDFX(l), PDFX(2), LMI, LPL, PDFX (-3) , 
PDFX(-4), PDFX(-5), PDFX(O) 

14 FORMAT (1PE6.0,1X,8(1PE11.4)) 

1 CONTINUE 
STOP 
END 
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Here is what this program returns. The main table has been shghtly edited (mainly 
E-07 — > E-7 etc.), to avoid having to use a yet smaller font. 

AUV = 0.333333 ADV = 0.137931 ALS = 0.136565 ASS = 0.027313 AGL = 0.364858 
NUV = 5.107200 NDV = 3.064320 NLS = 0.775950 NSS = 0.155190 NGL = 1.700000 

ALPHA_S(MR2(M2)) = 0.116032 
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The first line of the output provides the respective momentum fractions carried by the 
valence quarks, the light-quark sea and the gluons at the initial scale. This information 
has been obtained by accessing the common-block PANORM filled for this purpose by the 
inner input routine not described in this section. The second line are the corresponding 
normalization factors, cf. Eq. (4.1). The rest is a set of reference results which agree, 
except for some marginal offsets at x — 0.9 for the tiny sea-quark distributions, with the 
upper part of table 4 in ref. [37]. A test run by the user should lead to the same numbers. 

A final remark on the efficient calculation of the parton distributions at several scales. 
XPARTON saves the moments used for the inversion, and will re-use them the next time — 
recall that the a{N) in Eq. (3.2) do not know about x — unless the input or the scale 
//^ have been changed. Consequently one should first perform the Mellin inversion for all 
desired x-vahics at one scale before proceeding to the next scale, instead of ordering the 
calls of XPARTON in another manner. The same applies, of course, also to corresponding 
routines determining obscrvables like structure functions from the A^-space expressions. 

4.5 The longitudinally polarized case 

The evolution of the polarized parton distributions is performed completely analogous to 
the unpolarized case discussed above. Indeed, one could write the program such that the 
same user-interface routines would deal with both cases. Here we have decided to keep at 
least the polarized and unpolarized initialization and input routines separately. 

CALL INITPOL('EVOLPAR') 

initialises the polarizated evolution. The options are almost identical to those discussed in 
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section 4.1. The order of the evolution is restricted to NPORD = and 1 at present, since the 
polarized three-loop sphtting functions are not yet known. The respective files supplying 
the initialization parameter for EVOLPAR = 1 and 2 are usrpinit . dat and usrpinit . f 
providing the subroutine USRPINIT. The defaults are the same as in section 4.1. 

The input parameters and initial polarized distributions for the evolution are set by 

CALL INITPINP('INPPAR') . 

The corresponding data files are usrpinp . dat and usrpinp . f containing the routine 
USRPINP. Instead of Eq. (4.1) the default toy input for program checks reads 

xAu^ = +1.3 - x f {1 + 3x) 

xAd^ = -0.5 x^-^ (1 - x f (1 + 4x) 

xAg = +1.5 a;°-^ (1 - x f 

xAd = xAu = -0.05 (1 -a;)^ 

xAs = xAs = +0.5 xAd . (4.4) 

The other input parameter have to same meaning (and the same defaults) as in section 
4.2, with the exception of 

IMOMIN 

For IMOMIN = (used as the internal default here), we put A^j = 1 in Eqs. (3.7) and (3.8) 
for all seven input combinations u^, d^, L_ = d — u, L+ — 2{d + u), s± — s ± s and g. 
Otherwise all Ni are chosen such that the first elements of PUV, PDV, PLM, PLP, PSM, PSP 
and PGL represent the first moments of the respective initial distributions. PSM (5) is a 
normal input parameter in the polarized case. 



4.6 Output in AT-space 

The above setup assumes that the user is interested in the x-space parton distributions 
as obtained by the method outlined in section 3.1. For those who want to evolve fixed 
moments (e.g., momentum fractions) or prefer to invert the A^-space results by other 
means (e.g., Monte-Carlo integration), we also provide the user-interface subroutine 

NPARTON (PDFN, AS, 'N', 'M2', 'IPSTD\ 'IPOL', 'EVOLPAR', 'INPPAR') . 

There is no point to have separate overall and input initialization calls here, as all quanti- 
ties need to be recalculated each time anyway as N is a (double-complex) free parameter. 
Consequently the calls of the initialization and input routines have been integrated into 
NPARTON, which therefore takes over the respective parameters EVOLPAR and INPPAR (the 
options remain as before). In this case also the choice between the unpolarized (IPOL =0) 
and polarized evolution (otherwise) is made via the call of NPARTON. Note that this mode 
of using the program is an option for which the code has not really been optimized. The 
double-complex output array has to be declared as PDFN (-6: 6) in the calling program. 
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5 A short reference guide 



In this section we provide brief descriptions of the subroutines and functions included in 
the evolution package. Unless specified otherwise, the function or subroutine 'ROUTINE' 
is stored in the file 'routine.f, where more information can be found. This results in 
quite a few files, but it appears easier this way to keep track of future upgrades and user 
modifications. A very large part of the internal communications in the program proceeds 
via named common-blocks. A corresponding list is included at the end of this section. 

5.1 Initialization routines 

The main initialization routine for the A^-space based parton evolution is 

5.1.1 INITEVOLCIPARO 

The options for IPAR have already been discussed in section 4.1. This routine sets some 
constants like the Euler-Mascheroni constant 7e = ZETA(l), the lowest values of the 
Riemann's Zeta function, Ci>i = ZETA(i), and the lowest SU(3) invariants CA, CF and 
TR. It provides the fixed array NA of complex moments and the corresponding weights 
WN for the Gauss integrations discussed in section 3.1. The analytically continued simple 
harmonic sums Si{N) are calculated for i = 1, . . . , 6 on these support points and stored 
as S(K,i) where K is the parameter of the array NA. Finally the routine calls, depending 
on the initialization parameters discussed in section 4.1, the subroutines for the A^-space 
splitting functions, Uy-matching conditions and evolution matrices listed below. 

INITEVOL also sets a couple of internal initialization parameters. The consistency of 
the order NAORD of the coupling constant (sec section 2.1) with that of the parton evolution 
is enforced by NAORD = NPORD. The number of steps for the Runge-Kutta integration of 
Eq. (2.2) and the maximal power of Og in the [/-matrix solution (2.23) are set to 

NASTPS = 20 and NUORD = 15 . 

Note that NUORD is presently restricted by array declarations to values of 20 or less. 

5.1.2 USRINIT (IFAST, IVFNS, NFF, IMODEV, NPORD, FR2) 

The subroutine providing the initialization parameters if INITEVOL is called with IPAR= 2. 

5.1.3 BETAFCT 

Provides the values BETA(O), . . ., BETA (3) (2.3) of the QCD /3-function for = 3, . . . , 6. 

5.1.4 PSI(Z), DPSI(Z,M) in the file psifcts.f 

The complex ■0-function PSI(Z) and its m-th derivatives ■0^"*^-^) — DPSI(Z,M) calculated 
using the functional equations and the asymptotic expansions as discussed in section 3.2. 
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5.1.5 PNSOMOM, PSGOMOM in the file pnsgOmom.f 

The lowest-order non-singlet (routine PNSOMOM) and singlet (routine PSGOMOM) A'"-space 
splitting function PONS and POSG on the array NA with parameter KN defined by INITEVOL, 
as all splitting functions for = 3, . . . , 6. The singlet sector uses (also at higher orders) 
the matrix notation (2.11) with the array arguments 1 = q and 2 = g. 

5.1.6 PNSIMOM 

The subroutine for the NLO non-singlet splitting functions PINS at = on the array 
NA. The routine needs to be called before PSGIMOM, since it provides part of -Pqq\ see 
Eq. (2.15), and the non-trivial analytic continuations according to Appendix (A.l) of 
ref. [6] and Eq. (3.4). Starting from NLO the three cases (2.14) are included via an 
additional array dimension with arguments 1 = +, 2 = - and 3 = v (= - at NLO). 

5.1.7 PSGIMOM 

The corresponding subroutine for the NLO singlet splitting functions PlSG(KN,NF,i, j). 

5.1.8 PNS2M0M 

The NNLO non-singlet sphtting functions P2NS in A^-space, as always on the array NA, as 
obtained from the accurate a;-space parametrizations (4.22)-(4.24) of ref. [13] since the 
complex moments of the exact results are presently unknown. Analogous to the NLO 
case, this routine needs to be called before the singlet case PSG2M0M. Also here /^^ = A*- 

5.1.9 PSG2M0M 

The corresponding routine for the singlet quantities P2SG, using (4.32)-(4.35) of ref. [14]. 

5.1.10 LSGMOM 

The eigenvalue decomposition (2.26)-(2.28) of the LO singlet sphtting-function matrix 
POSG divided by BETAO. The A^- and n^-dependent eigenvalues are denoted by R(KN , NF , 1) , 
the corresponding projection operators by E(KN,NF,i, j ,1) with 1 = 1, 2. 

5.1.11 USGIMOM 

The subroutine for the NLO singlet evolution matrix Ui= U1(KN,NF,I,J) in A^-space 
obtained from Eq. (2.31) for k — 1. Here and in the following routines the terms arising 
from 7^ in Eq. (2.8) are included. This routine needs to be called before USGIHMOM. 
Procedures (presently absent) for scheme transformations of the splitting functions (pro- 
vided in MS by the above routines) would have to be called before this routine. 
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5.1.12 USGIHMOM 



Adds the higher-order contributions Uk= U1H(K, . . ) to the [/-matrices for the iterated 
NLO solutions of the evolution equations discussed in section 2.4. Note that there are 
no non-singlet routines corresponding to USGIMOM and USGIHMOM since Uf" are trivial and 
Eqs. (2.34) and (2.35) are used for the iterated non-singlet solutions instead of Eq. (2.23). 

5.1.13 UNS2M0M 

The routine for the NNLO non-singlet evolution operators UNS2, including the higher- 
order pieces for the iterative solutions. Also here the quantities for the evolution of q^"^ 
are included by an extra array dimension with the argument 1 = +, 2 = - and 3 = v. 

5.1.14 USG2M0M 

As the subroutine USGIMOM, but providing U2 for the truncated NNLO solution {k — 2). 

5.1.15 USG2HM0M 

As USGIHMOM, but adding the higher terms U2H required for the iterated NNLO solutions. 
The singlet [/-matrix routines partly build upon each other, hence they should be called 
in a proper order, like the one used in this list. 

5.1.16 ANS2M0M 

The (NNLO) non-singlet coefficient A^^^f^ = A2NS(KN) for the MS flavour-number 
transition (2.39) at //^ — m\, obtained from the x-space results in Appendix B of ref. [21]. 

5.1.17 ASG2M0M 

The flg (NNLO) A^-spacc singlet coefficients A2SG in Eq. (2.40) obtained from the same 
source. The routine uses A2NS, hence it is to be called after ANS2M0M. The moments of 
the parametrization (3.5) are used for A^^^'^^ A2SG (KN , 1 , 2) . 

5.2 Input and flavour- threshold routines 

The steering routine for the iV-space distributions at initial scale and flavour thresholds is 
5.2.1 INITINP('IPAR') 

The options switched by IPAR have already been discussed in section 4.2. This routine calls 
the appropriate routine INPLMOMl or INPLM0M2 and calculates the initial value of for the 
evolution, which is unequal to the (properly normalized) input parameter in section 4.2 
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for 7^ n. For IVFNS ^ 0, EVNFTHR is then used to store also the parton distributions 
at the heavy-flavour thresholds //^ = m^. This is another efficiency measure: Recall that 
the evolution is performed in fixed-n^ steps, using Eqs. (2.39)-(2.41) in between. For all 
values /i^ > ml, e.g., the evolution to /i^ — ml is thus the same for a given input, and 
there is no point to repeat this part of the computation for every new value of //^. 

5.2.2 USRINP (PUV, . . . ,PGL, M20, ASI, MC2,MB2,MT2, NFORM, IMOMIN, ISSIMP) 
The subroutine providing the input parameters if INITINP is called with IPAR= 2. 

5.2.3 INPLMOMl (PUV, PDV, PLM, PLP, PSM, PSP, PGL, IMOMIN, ISSIMP) 

This routine calculates, from the input parameters discussed already in section 4.2, the 
moments of the ansatz (3.7) on the array NA and stores them in the basis (2.16) with the 
notation ql^{N) = VAI(KN), = SGI, = M'k'I, -y^f = P'k'I and ^ = GLI for A; = 3, 8. 
Also the common-block with the second moments and normalizations shown in section 
4.4 is written here. The routine sets the flag IINNEW = 1 to inform other routines, like 
XPARTON discussed in 4.3, that a new input call has taken place. 

5.2.4 INPLM0M2 (PUV, PDV, PLM, PLP, PSM, PSP, PGL, IMOMIN, ISSIMP) 

As the previous routine, but for the ansatz (3.8) for the three-flavour initial distributions. 
Other input forms can be added as INPLM0M3 etc. if needed, which requires only a miminal 
additional modiflcation of INITINP in section 5.2.1. 

5.2.5 EBETA(Z1,Z2) in the file ebetafct.f 

The complex Beta function EBETA(Z1,Z2) calculated from the asymptotic expansion of 
lnr(2;) as discussed in section 3.2. Called (only) by the routines for the A^-space inputs. 

5.2.6 EVNFTHR (MC2, MB2, MT2) 

For IVFNS = 1, this routine is used to evolve Og and the partons from the three-flavour 
initial scale to the four- to six-flavour thresholds /^^ = m^. The results for are stored 
as ASC, ASB and AST, those for ql^ as VAC, VAB and VAT, etc. For ml > 10^° or mg^ > 10^° 
(in GeV^, as always), the corresponding parts of the calculations are skipped. 

5.2.7 ASNFl (ASNF, LOGRH, NF) in the file asmatch.f 

This functions computes as ^'^ ^ = ASNFl from as*' = ASNF and ln(/i^/m^) = LOGRH 
according to Eqs. (2.41)-(2.42) at maximally N^LO. The order of the n^-matching is set 
by NAORD specifled via the initiahzation parameter NPORD, see section 5.1.1. 
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5.3 Evolution and Mellin-inversion routines 

5.3.1 AS (R2, R20, ASO, NF) in the file asrgkt.f 

Returns Og = AS at the renormalization scale /i^ = R2 as obtained beyond LO from a 
fourth order Runge-Kutta integration of Eq. (2.2) in NASTPS steps from the initial value 

00 = ASO at /x^o = ^^20 at order NAORD for NF quark flavours. At LO Eq. (2.4) is used. 
Possible other versions have to be provided (in the same notation) by the use her- /himself. 

5.3.2 ENSGON (ENS, ESG, ASI, ASF, S, KN, NF) 

The subroutine for the LO non-singlet and singlet iV-space evolution operators (2.22) at 
a fixed number of flavours NF. The respective kernels ENS and ESG(I, J), I, J = 1, 2 with 

1 = q, 2 = g, are returned for a moment N specified via the counter KN on the array 
NA. The initial and final scales arc specified by the respective values Oq = ASI and = 
ASF of the coupling constant. S = ln(ao/as) is not calculated internally for efficiency. 

5.3.3 ENSIN (ENS, ALPI, ALPF, S, KN, NF) 

This routine returns one of the NLO hadronic non-singlet kernels ENS(K) (2.33)-(2.36) 
for the evolution of the '+' (K=l) and — = v (K=2,3) combinations (2.13) of quark 
distributions. The mode of the evolution is set via the initialization parameter IMODEV. 
The input arguments of the routine are as described for the previous routine. 

5.3.4 ESGIN (ESG, ASI, ASF, S, KN, NF) 

As ENSIN, but for the singlet kernels ESG (I, J) according to Eqs. (2.23) and (2.24). 

5.3.5 ENS2N (ENS, ASI, ASF, S, KN, NF, NSMIN, NSMAX)) 

As ENSIN, but for the NNLO kernels ENS(K) using Eqs. (2.23), (2.24) and (2.33). The 
additional arguments NSMIN and NSMAX set the range in K for which ENS (K) is defined. This 
allows to save time, e.g., in e.m. structure- function calculations requiring only ENS(l). 

5.3.6 ESG2N (ESG, ASI, ASF, S, KN, NF) 

As ESGIN, but for the NNLO kernels ESG (I, J) governing the evolution of the flavour- 
singlet parton distributions. Recall that all these kernels assume a fixed ratio /i/jir- 

5.3.7 EVNVFN (PDFN, ASI, ASF, NF, NLOW, NHIGH, IPSTD) 

The subroutine providing the TV-space parton distributions PDFN for the variable flavour- 
number evolution. These results are returned for the part NLOW < KN < NHIGH of the 
array NA in a notation chosen via IPSTD as discussed in section 4.3. The routine makes 
use of the flavour-threshold results stored by EVNFTHR described in section 5.2.6. 
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5.3.8 EVNFFN (PDFN, ASI, ASF, NF, NLOW, NHIGH, IPSTD) 



As EVNVFN, but for the fixed flavour-number evolution with 3 < NF < 5 partonic flavours. 

5.3.9 XPARTON (PDFX, AS, 'X^ 'M2', ' IFLOW , 'IFHIGH', 'IPSTD') 

The top-level evolution and Mellin-inversion routine for standard use of the program, 
discussed already in section 4.3. It checks whether the previous results for Og and PDFN can 
be re-used (same M2 and no new input call), calculates the otherwise required quantities, 
calls EVNVFN or EVNFFN if necessary, and carries out the Gauss integrations. 

5.4 Routines for the polarized case 

To avoid a substantial duplication of code at this point, the polarized evolution uses the 
internal routines described above as far as possible. Thus a simultaneous evolution of 
unpolarizcd and polarized parton distributions, e.g., for A^-space analyses ofpp scattering 
with only one proton polarized, is not possible at this point. The additional routines are: 

5.4.1 INITPOL('IPAR') 

The polarized version of the main initialization routine INITEVOL in section 5.1.1. The 
same options are available except for the restriction of the evolution to LO and NLO. 

5.4.2 USRPINIT (IFAST, IVFNS, NFF, IMODEV, NPORD, FR2) 

The subroutine providing the input parameters for INITPOL if called with IPAR= 2. 

5.4.3 PSGOPMOM 

The lowest-order polarized singlet A^-space splitting function POSG corresponding to the 
unpolarized routine in section 5.1.5. Note that there is no corresponding non-singlet 
routine, as in this case the unpolarized and polarized splitting functions are identical. 

5.4.4 PNSIPMOM, PSGIPMOM 

The polarized counterparts of PNSIMOM and PSGIMOM described in sections 5.1.6 and 5.1.7. 
The same array names and arguments are employed for use by the general [/-matrix and 
evolution routines in sections 5.1 and 5.3. The non-singlet quantities are identical to the 
unpolarized case, but with the and - = v entries interchanged. 

5.4.5 ANS2PM0M, ASG2PM0M 

The subroutines for the polarized case corresponding to ANS2M0M and ASG2M0M in sections 
5.1.16 and 5.1.17. Presently these are dummy routines. The presence of the arrays A2NS 
and A2SG is technically required in EVNFTHR (section 5.2.6). 
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5.4.6 INITPINPCIPARO 



The input initialization routine for the polarized evolution corresponding to INITINP in 
section 5.2.1. I PAR switches the input options as discussed in section 4.2 and 4.5. 

5.4.7 USRPINP (PUV, . . . ,PGL, M20, ASI, MC2,MB2,MT2, NFORM, IMOMIN, ISSIMP) 
The subroutine providing the input parameters for INITPINP if called with IPAR= 2. 

5.4.8 INPPMOMl, INPPM0M2 

The input- moment subroutines corresponding to INPLMOMl and INPLM0M2 in sections 5.2.3 
and 5.2.4. The only difference of the input options has been discussed in section 4.5. The 
common-block PANORM now returns the first instead of the second moments as AUV etc. 

5.5 Routines for output in AT-space 

The user interface for obtaining the results of the evolution at any complex value of N is 

5.5.1 NPARTON (PDFN, AS, 'N' , M2', 'IPSTD', 'IPOL', 'EVOLPAR', 'INPPAR') 

briefly discussed in section 4.6. PDFN (-6: 6) is the double complex output array with two 
options switched by IPSTD as explained for the analogous x-space array PDFX in section 
4.3. AS returns the corresponding value as(/x^(/i^)) with /x^ = M2 (both double-precision). 
Note that NPARTON is called without previous initialization and input calls as these calls 
are issues by this routine using the values of IPOL, EVOLPAR and INPPAR. 

5.5.2 INITMOM (IPOL, EVOLPAR) 

The special initialization routine called by NPARTON. Depending on IPOL this routine takes 
over the respective roles of INITEVOL in section 6.1.1 or INITPOL in section 6.4.1 in the 
calculations for one moment A^. The array dimension NA for the standard setup is formally 
kept in order to comply with the fixed declarations in the low-level routines. 

5.6 A list of common blocks 

Here we present the list of common blocks mentioned above, ordered alphabetically. The 
dimensions of the array variables are not indicated for brevity. Most, but not all of the 
variables stored by these common blocks have been mentioned above. 

COMMON / ANS2 / A2NS 

COMMON / ASFTHR / ASC, M2C, ASB, M2B, AST, M2T 
COMMON / ASG2 / A2SG 



29 



COMMON / ASINP / ASO, M20 

COMMON / ASPAR / NAORD, NASTPS 

COMMON / BETA / BETAG, BETAl, BETA2, BETAS 

COMMON / COLOUR / CF, CA, TR 

COMMON / EVMOD / IMODEV 

COMMON / FRRAT / LOGFR 

COMMON / HSUMS / S 

COMMON / INPNEW / IINNEW 

COMMON / INVEST / IFAST 

COMMON / ITORD / NUORD 

COMMON / KR0N2D / D 

COMMON / LSG / R 

COMMON / MOMS / NA 

COMMON / NCONT / C, CC 

COMMON / NFFIX / NFF 

COMMON / NFUSED / NFLOW, NFHIGH 

COMMON / NNUSED / NMAX 

COMMON / ORDER / NPORD 

COMMON / PABTHR / VAB, MSB, MSB, M15B, M24B, SGB, PSB, P8B, P15B, 

P24B, GLB 

COMMON / PACTHR / VAC, MSC, M8C, M15C, SGC, PSC, P8C, P15C, GLC 
COMMON / PAINP / VAI, MSI, M8I, SGI, PSI, P8I, GLI 

COMMON / PANORM / AUV, ADV, ALS, ASS, AGL, NUV, NDV, NLS, NSS, NGL 
COMMON / PATTHR / VAT, M3T, MBT, M15T, M24T, MS5T, SGT, PST, P8T, 



P15T, P24T, PS5T, GLT 
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/ 
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6 Accuracy and speed of the program 



The main numerical step in an iV-space evolution program is the inverse Mellin trans- 
formation of the final results back to x-space. As discussed in section 3.1, the present 
program uses a one-fits-all grid of fixed A^-values for optimal performance in large com- 
putations. The accuracy of the programmed inversion can be easiest tested by calling 
XPARTON at this initial scale M2 = /Iq and dividing by the known x-space initial distri- 
butions. This has been done in Fig. 2 for the two extreme shapes in Eq. (4.1), and 
= u. Recall that xu^ vanishes as x^'^ for x — *• 0, while xu^ ~ (1 — xV for x — > 1. 
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Figure 2: The up-valence and up-sea distributions at the initial scale (4.2) obtained by 
the standard and fast Mellin inversions, divided by the exact x-space inputs in Eq. (4.1). 
The same four curves are shown in the left (right) part for small (large) values of x. 



The results for the standard inversion using up to 144 moments (IFAST = 0) for this 
test do not show any noticeable deviations from the exact values at 10~^ < x < 0.9 for m^, 
and at 10~^^ < x < 0.8 for u^- These ranges refer to the program's default value c = 1.9 
of the contour abscissa in Eq. (3.2). Smaller values would improve the inversion of at 
extremely small x at the cost of worsening that of Us at very large x. 

The faster, but close to the end-points inevitably less reliable option using between 32 
and 80 moments (IFAST = 1) can be safely employed at least for 10""^ $ x < 0.7, where 
the upper limit includes the effect of the further softening of the sea-quark distributions by 
the evolution. In many applications the small (and anyway poorly known) distributions 
XMv at smaller x and Ms at larger x are practically irrelevant; then the fast inversion can 
be used over the full range of x shown in the figure. 
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At NNLO the present overall accuracy of the program at 'normal' values of x is 
dominated by that of the parametrizations of the corresponding splitting functions and 
the riy^-matching coefficient discussed in section 3.2. An update of the NNLO benchmark 
results in ref. [37] with the complete three-loop splitting functions of refs. [13, 14] will be 

presented elsewhere [39]. For checks of tables 5 and 6 in ref. [37] the user needs to employ 
the previous approximations [36] instead of the routines 5.1.8 and 5.1.9. For the moment 
these approximations are still available as pns2mold . f and psg2mold . f . 

Finally we illustrate the speed of the program. For this purpose the input (4.1) is 
evolved, using Eqs. (4.2) and (4.3), for variable to 20 scales ^ and the results are 
Mellin-inverted at 25 values of x at each scale. The chosen scales and values of x are 



DATA MS 


/ 


2 


.ODO, 


2.7D0, 


3.6D0, 


5. DO, 


7. DO, 


l.Dl, 


1 




1, 


.4D1, 


2.D1, 


3.D1, 


5.D1, 


7.D1, 


1.D2, 


2 




2 


.D2, 


5.D2, 


1.D3, 


3.D3, 


1.D4, 


4.D4, 


3 




2 


.D5, 


1.D6 / 










DATA XB 


/ 


1, 


.D-8, 


l.D-7, 


l.D-6, 








1 




1 


.D-5, 


2. D-5, 


5. D-5, 


l.D-4, 


2.D-4, 


5.D-4, 


1 




1 


.D-3, 


2. D-3, 


5. D-3, 


l.D-2, 


2.D-2, 


5.D-2, 


2 




1, 


.D-1, 


1.5D-1, 


2. D-1, 


3. D-1, 


4. D-1, 


5. D-1, 


3 




6, 


.D-1, 


7. D-1, 


8. D-1, 


9. D-1 / 







The whole calculation, including the input initialization but excluding the general ini- 
tiahzation, is then repeated 200 times. Thus we perform precisely 10^ calls of XPARTON. 
The times (in seconds) required for this computation on a 2.0 GHz Pentium-IV processor 
are listed below for various values of NPORD, IMODEV and IFAST (defined in section 4.1). 
The value of the latter parameter is given after the name of the compiler. We have used 
the GNU compiler g77 (version 3.2-36) with -02 and the (non-commercial) 'Intel Fortran 
Compiler 8.0 for Linux' if ort with -xW (optimization with -02 is the default for if ort). 



Type of evolution 


g77,0 


g77,l 


if ort, 


if ort, 1 


LO 


7.8 


4.3 


2.5 


1.4 


IMODEV=0, NLO 


8.2 


4.5 


2.7 


1.5 


IMODEV=0, NNLO 


9.2 


5.0 


3.2 


1.8 


IM0DEV=1, NLO 


10.7 


5.8 


3.8 


2.1 


IM0DEV=1, NNLO 


11.3 


6.2 


4.4 


2.4 



Thus roughly 10^ to 10^ Melhn inversions of the complete set of partons can be per- 
formed per second under the above conditions, depending on the initialization parameters 
and the available compiler. As expected, the evolution becomes somewhat slower with in- 
creasing order, and the iterative solutions (involving high orders of the evolution matrices 
[/) are slower than the truncated evolution. 
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7 Summary and outlook 



We have presented the fast, flexible and accurate Fortran package QCD-Pegasus 
solving the evolution equations for the unpolarized and polarized parton distributions 
of hadrons in perturbative QCD. The code is designed such that the program is relatively 
fastest when this is most needed, viz in large repetitive computations like fits of the initial 
distributions or error analyses according to refs. [40, 41, 42]. In standard applications, 
the user needs to interact with only very few top-level routines which communicate the 
evolution parameters and initial distributions to, and the results of the evolution from 
the lower-level routines actually performing the underlying evolution in Mellin-A^ space. 

The full potential of the A^-space approach is realized once not only the parton evolu- 
tion, but also the convolutions with the partonic cross sections (coefficient functions) are 
performed (as multiplications of the moments) before Mellin-inverting back to x-space. 
Doing this poses no problems if the coefficient functions are known (to a sufficient accu- 
racy) in a form facilitating the calculation of their complex-A^ moments. Even if this is 
not the case — for example, if coefficient functions are only known via an x-space Monte- 
Carlo programs including experimental cuts — there are efficient methods to proceed in 
iV-space, see ref. [43] and especially refs. [44, 45]. Sample routines for the latter proce- 
dure will be presented elsewhere for representative observables in ep and pp/pp scattering. 
The user can then implement the processes he or she is interested in along these lines. 
QCD-Pegasus can of course also be used efficiently within 'normal' x-space analyses. 

If additional routines are included to deal with the additional 'pointlike' solution, the 
program can also be used for the QCD evolution of the parton distributions of the photon. 
These routines largely exist, but have not been included in this first public version of the 
program. An even easier extension would be the evolution of also 'timelike' distributions 
(fragmentation functions) of effectively massless partons. Considerably larger modifica- 
tions would be required to include other (hypothetical) light coloured particles like a light 
gluino which would lead to a enlarged 3x3 flavour-singlet sector. This extension was 
actually present in an earlier incarnation of this program, see ref. [46], but has not been 
pursued anymore since the light-gluino window is now generally considered closed. 

The code of QCD-Pegasus can be obtained by downloading the source of this article 
from the hep-ph preprint archive. Relevant changes of the code, if required, will be com- 
municated by replacing the article in this archive. The program is distributed under the 
GNU public license, with the obvious additional requirement that scientific publications 
(including preprints) which have been prepared using the present code, or a modified 
version of it, should include a reference to the present article. The author can be reached 
by electronic mail for bug reports, comments and questions. 
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